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ABSTRACT 



5_i 

We present an idealized, spherical model of the evolution of a magnetized molecular 
cloud due to ambipolar diffusion. This model allows us to follow the quasi-static 
evolution of the cloud's core prior to collapse and the subsequent evolution of the 
remaining envelope. By neglecting the thermal pressure gradients in comparison with 
magnetic stresses and by assuming that the ion velocity is small compared with the 
neutral velocity, we are able to find exact analytic solutions to the MHD equations. We 
^-j. ; show that, in the case of a centrally condensed cloud, a core of finite mass collapses into 

the origin leaving behind a quasi-static envelope, whereas initially homogeneous clouds 
never develop any structure in the absence of thermal stresses, and collapse as a whole. 
Prior to the collapse of the core, the cloud's evolution is characterized by two phases: 
a long, quasi-static phase where the relevant timescale is the ambipolar diffusion 
time (treated in this paper), and a short, dynamical phase where the characteristic 
timescale is the free-fall time. The collapse of the core is an "outside-in" collapse. The 



O 

quasi-static evolution terminates when the cloud becomes magnetically supercritical; 
thereafter its evolution is dynamical, and a singularity develops at the origin — a 
protostar. After the initial formation of the protostar, the outer envelope continues 
to evolve quasi-statically, while the region of dynamical infall grows with time — an 
"inside-out" collapse. We use our solution to estimate the magnetic flux trapped in 
the collapsing core and the mass accretion rate onto the newly formed protostar. 

Our results agree, within factors of order unity, with the numerical results of Fiedler 
& Mouschovias (1992) for the physical quantities in the midplane of a collapsing, 
magnetized, axisymmetric cloud up to the onset of dynamical collapse. Our simple 
approach thus captures the basic physics of a self-gravitating, magnetized cloud in 
which the evolution is driven by ambipolar diffusion. It also enables us to treat 
the evolution of the accretion onto the protostar after collapse, for which detailed 
numerical results are as yet unavailable. Remarkably, we find that at late times the 
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accretion rate becomes comparable to that of a non-magnetized, singular isothermal 
sphere, provided that the ionization is due to galactic cosmic rays. 

Subject headings: diffusion — hydromagnetics — ISM: clouds — ISM: magnetic fields — 
plasmas — stars: formation 

1. Introduction 

How do stars form? This age-old question is one of the outstanding challenges of modern 
astrophysics. Even though we still lack a definite answer to this question, some partial answers 
have begun to emerge over the last three decades. The birth of molecular radio astronomy and 
infrared astronomy enabled us to definitively identify molecular clouds as the site of star formation, 
and to probe the physical conditions therein (Zuckerman & Palmer, 1974). Notwithstanding the 
spectacular advances in instrumentation of the last decade and the resulting explosive growth in 
the number of infrared and radio observations, the detection of a cloud undergoing dynamical 
collapse is very difficult; to date there are only a few, tentative such detections (see Zhou et al., 
1994 and references therein). By necessity, then, one has to resort to theory to understand how 
the molecular gas ultimately condenses into a star. 

The efforts in this direction can be traced back to the pioneering work of Larson (1969) and 
Penston (1969). In these early calculations only thermal pressure — and sometimes rotation — were 
included, while magnetic fields were neglected (see Mouschovias [1978] and Bodenheimer k, Black 
[1978] for a list of early references). Of particular importance are the semi-analytic similarity 
solutions of Shu (1977) and Hunter (1977). These spherical solutions, though idealized, shed 
important light on the nature of the non-magnetic collapse. Non-magnetic collapse calculations 
are still pursued today; the current state of development in this field is well represented by the 
work of Myhill & Boss (1993), Boss & Myhill (1993), and Foster & Chevalier (1993). 

The relevance of magnetic fields to star formation were first pointed out by Mestel & Spitzer 
(1956). Using the virial theorem, Strittmater (1966) calculated the critical mass for gravitational 
collapse perpendicular to the field lines of a cold, uniformly magnetized cloud, and Field (1965) 
studied the effect of the magnetic field on the thermal instability. These studies established that 
the presence of a magnetic field introduces qualitative changes to the conclusions of non-magnetic 
studies. However, not until the late seventies were detailed numerical calculations of magnetic 
clouds carried out. Mouschovias (1976a, b) performed the first self-consistent numerical studies of 
equilibrium, self-gravitating magnetic clouds, and Scott & Black (1978) were the first to study the 
dynamics of collapsing, non-rotating clouds with a frozen-in magnetic field. With the introduction 
of a magnetic field, a realistic description requires at least two spatial dimensions. Moreover, at 
the densities of interest (n ^ 10 3 cm -3 ) the degree of ionization is low, and one has to consider 
the ion and neutral material as two separate fluids. This two-fluid nature of the problem has 
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important consequences for the stability and collapse of magnetized clouds. The drift of the 
neutrals relative to the ions — a process known as ambipolar diffusion (Mestel & Spitzer, 1956) — is 
responsible for bringing to the verge of collapse an otherwise stable cloud. As ambipolar diffusion 
proceeds, the density of the central region grows until the inner core attains a mass larger than a 
critical mass. Once this happens, the central region of the cloud is in a state of dynamical collapse 
(for a review, see Mouschovias [1996] and McKee et. al [1993] and references therein). 

Except for some idealized, one-dimensional solutions which illustrated the role played by 
ambipolar diffusion in the redistribution of magnetic flux (Mouschovias 1979; Shu 1983; Scott, 
1984), the study of the stability and collapse of realistic, magnetized clouds has been the domain 
of complex numerical calculations, with major contributions by Nakano (1979, 1982, 1983); 
Mouschovias and collaborators (Mouschovias, Paleologou, & Fiedler 1985; Mouschovias & Morton, 
1991; Fiedler & Mouschovias, 1992; Basu & Mouschovias, 1994, 1995a,b; Ciolek & Mouschovias, 
1993, 1994a,b, 1995); Tomisaka, Ikeuchi, & Nakamura (1988a,b; 1989; 1990); Tomisaka, (1991); 
Lizano & Shu (1989); and Galli & Shu (1993a,b). In contrast to the non-magnetic models, none 
of these calculations can follow the cloud's evolution up to and beyond the point where a star is 
formed. They either stop at an early stage of dynamical collapse (Fiedler & Mouschovias, 1992); 
or tackle the problem after the formation of the star and follow the infall of the envelope under 
the gravitational field of a point source (Galli & Shu, 1993b). 

The shortcomings of current numerical work, and the lack of analytic or semi-analytic models 
to complement and illuminate the numerical results, motivated us to undertake the present study. 
Our goal is to understand — by means of a very simple model — how a globally stable cloud evolves, 
through ambipolar diffusion, into a state where a central core collapses, and to shed some light 
on the evolution of the remnant envelope. By doing so, we hope to distill the basic features that 
can be lost in the detail and complexity of numerical calculations. Previous attempts along these 
lines adopted either planar (Shu, 1983; Paleologou & Mouschovias, 1983) or cylindrical symmetry 
(Mouschovias & Morton, 1991). These clouds are always magnetically subcritical under conditions 
of flux freezing, and are of limited relevance to the study of collapsing clouds (Mouschovias, 1991; 
McKee et al., 1993). In our case, we adopt spherical symmetry, where the only non- vanishing 
component of the Lorentz force is the scalar term oc dB 2 /dr. In contrast to the planar and 
cylindrical cases, we cannot identify the actual magnetic field lines, but a three-dimensional cloud 
can undergo collapse even when flux freezing holds (Mouschovias, 1991; McKee et al., 1993). 

In this paper, we neglect thermal pressure altogether, and set the gas temperature equal to 
zero. The case where thermal support is included is deferred to a second paper. Our motivation 
for this drastic approximation is twofold. From a practical viewpoint, this simplification enables 
us to treat ambipolar diffusion from the onset of cloud contraction to the collapse of a stellar core 
and beyond, an intrinsically non self-similar problem. Second, we show explicitly (§5.2.2) that 
the role of thermal pressure diminishes in time as the central density rises. Thus, our calculation 
should give an increasingly accurate representation of the cloud's later evolution, as it forms a 
protostar and then slowly drains onto the protostar. 
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The paper is organized as follows: in §2 we introduce our basic equations and approximations 
and present a general solution; the general properties of the solution for a particular family of 
initial conditions are presented in §3; in §4 we study the collapse of a homogeneous cloud; §5 
contains our results for a centrally condensed cloud and a comparison between our results and the 
numerical results of Fiedler & Mouschovias (1992); and our conclusions are presented in §6. 



2. Formulation and a General Solution 

We start with the general, two-fluid MHD equations governing ambipolar diffusion in 
self-gravitating molecular clouds (see, e.g., Paleologou & Mouschovias, 1983). We assume 
steady-state ionization — a situation that obtains in dense clouds — and adopt the ionization law 

n i = K i n l r / 2 , (1) 

where rij and n n are, respectively, the ion and neutral number density sand K{ is a constant with 
value Ki = 9.5 x 10 -6 cm" 3 / 2 when HCO + is the most abundant ion (Elmegreen, 1979). 

Our first basic approximation is that thermal stresses are negligible compared to magnetic 
forces, and we set the temperature of the gas equal to zero. The physical regime where this 
approximation is likely to be appropriate can be evaluated with the help of the parameter 

^Svr^T (2) 

where k-g is Boltzmann's constant, T is the gas temperature, and B is the strength of the magnetic 
field. For typical conditions in a starless cloud (see, e.g., Cernicharo, 1991; Heiles, 1987; Myers, 
1985) 

Therefore, low-mass clumps are magnetically dominated prior to the point at which their density 
is high enough for them to be visible as ammonia cores. The line widths of more massive clumps 
are often supersonic, so thermal pressure is relatively small there as well. However, we have not 
explicitly included turbulent pressure in our models because, in this formulation, its effects could 
be incorporated only under the artificial assumption that it is proportional to the static magnetic 
pressure. 

To further simplify the problem we impose spherical symmetry. The accuracy of this 
approximation is unclear. Theoretical models of cloud cores suggest that the cores are flattened 
along field lines (see, eg., Fiedler & Mouschovias, 1992; Lizano & Shu, 1989; Tomisaka, Ikeuchi, & 
Nakamura, 1988a). The observational evidence on cloud shapes has been interpreted as indicating 
that they are prolate (Myers et. al, 1991; Ryden, 1996) or that they are indeed flattened along 
field lines but that they are toroidal (Li & Shu 1996). Physically, the thermal and turbulent 
pressure, which we have neglected, prevent actual cloud cores from becoming too flattened, and 



- 5 - 




therefore we expect that our model will be qualitatively correct. Spherically symmetric models 
of magnetized clouds have been shown to represent the conditions for the onset of gravitational 
collapse reasonably well (Mouschovias & Spitzer, 1976; Holliman & McKee, 1996). Because we 
impose spherical symmetry, we cannot identify the actual magnetic field lines, nor can we account 
for nonradial forces. However, since we shall follow only the quasi-static part of the collapse 
(see below), we need the force equation only to estimate the flux remaining in the cloud, not to 
determine the dynamics of the collapse. As a result, we can adopt a simple approximation for the 
magnetic force, keeping the pressure gradient associated with the tangential fields but ignoring 
the associated tension. Let = Bq + B^ be the tangential field. The radial component of the 
magnetic force is then 

The last term represents the hoop stress, an inward force due to magnetic tension. However, 
in an open field configuration, which we assume, this force is more than compensated by the 
second term, which represents an outward tension. Since this second term cannot be estimated 
within the context of a spherical model, we simply ignore the last two terms, representing tension. 
Numerical simulations (cf. Fiedler & Mouschovias, 1992; Lizano & Shu, 1989) show that during 
the quasistatic phase the field does not develop a strong curvature. Therefore, our neglect of 
tension should be accurate within a factor ~ 2. 

To repeat, then, we ignore the last two terms in equation (Q) and keep the first term, 
representing pressure. We eliminate the angular dependence of B 2 by using its angular average 
(cf. Chiueh & Chou 1994); and, for simplicity, we shall use the notation B instead of {B 2 ) 1 / 2 , 
where ( ) stand for the angular average. 

Finally, at the typical densities of dense molecular clouds the ionization fraction is very small, 
and therefore, to a good degree of approximation, we have that 

dvi dvA 

_■> ~ 0. (6) 



and 



With these assumptions of negligible thermal pressure and of spherical symmetry, and 
neglecting the inertia of the ions, the momentum equations for the neutral and ion fluids and the 
mass conservation equation for the neutrals take the form 

Dv n {av) ni GMp. 

Dt m,i r 



rrij o7r or 



and 

Dpn Pn 9 (r 2 v n ) 
Dt r 2 dr 



0) 
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where 

D d d 

Dt = ai + V »* (10) 

indicates a substantial derivative; vd = Vi — v n is the ion-neutral drift velocity; {crv) ni is the 
average collision rate between neutrals and ions, with the value 1.7 x 10 _9 cm~ 3 s for HCO + -H2 
collisons (McDaniel & Mason, 1973); and all the other symbols have their usual meanings, with 
the subscripts i and n standing, respectively, for ion and neutral quantities. In writing equations 
@j and (jH), we have neglected the neutral mass in comparison with the ion mass, m n <C mi. We 
define a new constant k n i by 

k m ee = 8.2 x 10- 3 cm 3 / 2 s- 1 g- 1 / 2 , (11) 

where m n and mn are, respectively, the neutral's and proton's mass, and where we have set 
m n = 2.33 mn. 

At this point it is important to make contact with previous work by other authors. We choose 
to focus on the solutions by Mouschovias and collaborators (see references in §1) because our 
approach is closest to theirs. Their solutions depend on three dimensionless parameters, one of 
which is (Fiedler & Mouschovias, 1992) 

U& V37T 2 ) t ni fl ^ 12 ^ 

where 

1/2 



t '-'=(^) (13> 

is the free-fall time for their initial reference state with uniform neutral density p n q, and where 

tnifl = 7^-^ (14) 

is the neutral-ion collision time for the initial reference state (note that depends only on 
universal constants and the microphysics of the problem, and that we have again assumed 
m n <C mi). Our parameter k n i is directly related to us through 

k ni = 3.54G 1 / 2 z> fr . (15) 

Using k n i and equation (|T|), the system of eqs. (0) and (||) take the form 

Dv n , 3/2 GMp n 
Pnj^— = KiPrl V D — , (16) 

= - kniP Tv D -l^ (17) 



In principle, one has to add to eqs. (16), and (|T^) the induction equation to solve 
self-consistently for the magnetic field. We postpone our discussion of the induction equation until 
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we consider other approximations that apply to the solution of eqs. @, (jig), and (|17[) ; it turns 
out that, to the same order of approximation for which our solution of eqs. (0), (|i~6|), and ( |l7D 
holds, the induction equation is automatically satisfied to sufficient accuracy. 

Our fundamental approximations to solve eqs. (^), ([R]), and ( |17| ) are two. First, we assume 
that the evolution is quasi-static, and set 

5^-0. (18) 

Second, we assume that the ions are effectively static, so that their velocity is small compared 
even with that of the neutrals: 

Vi < v n , (19) 

so that 

vd ~ -v n - (20) 

These approxmations require comment. Recall that we are interested in the case where only 
the magnetic field supports the cloud against gravity. Because only the ions are directly coupled to 
the field, the neutrals are subject to the magnetic forces only through collisions with the ions, and 
there is a net slip of the neutral material past the magnetic field lines (Mestel & Spitzer, 1956). 
If this collisional coupling is effective, the drift velocity vjj will be small. There are two ways in 
which this can occur: On the one hand, when m ~ v n , the drift velocity is small, but the ions are 
essentially comoving with the neutrals; in other words, the cloud is not magnetically supported 
and it undergoes dynamical collapse (Mestel & Spitzer, 1956; note that when Vi = v n the first 
term on the rhs of equation ]l(| vanishes, and the material is in free-fall). On the other hand, vd 
is small also when < v„, and v n itself is small. In this case, the ions — and the magnetic field 
lines — are essentially stationary, and the neutral material slowly drifts by them. Thus, provided 
that v n is much smaller than the local free-fall velocity, the coupling between the ions and the 
neutrals is good, and the magnetic field is effective in supporting the cloud on timescales of order 
the free-fall time, but on a longer timescale the cloud is in a state of quasistatic contraction (e.g., 
Lizano & Shu 1989; Mouschovias, 1991). This is the case we consider here: quasi-static contraction 
of the neutrals through essentially static ions. 



Using eqs. ( |l~8| ) and (p0|), equation (16) becomes 

GM 



7 2 1/2- 
Kni f Pn 

Substitution into equation (|9|) then gives 

Dp n G p n d ( M 



Dt ~ k ni r 2 dr \ p VzJ ' 



(21) 



(22) 



Lagrangian coordinates are better suited for the solution of equation (|22|). By using the 
transformation 

— = (23) 
dM 4vrr 2 pn 
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on equation (22), we finally get our fundamental equation: 



Note that equation (p4|) is a self-contained equation for p n (t,M); once a solution is found, v n 
and r can be found from eqs. (|2l| ) and (^), respectively, while B is obtained by quadrature of 

which follows from eqs. (|l6|) and ( [T7D in the limit where eqs. (||), (||), and ( |T8| ) are valid. 



2.1. Self-Consistency and the Induction Equation 

One can check, a posteriori, the validity of the quasi-static approximation (eq. |l8|]) in the 
following way. For a given solution, evaluate the Lagrangian time derivative of v n . Then, in terms 
of the parameter a defined by 

Dv n _ GM 

DT = - a 7WY (26) 

the quasi-static approximation is valid as long as a <C 1. 

Our second approximation in solving eqs. (|9|),(|l6|), and (O) is given by equation (|l9|), i.e., 
that the motion of the ions is negligible compared to that of the neutral material. The magnitude 
of V{ is critical to the solution of the induction equation, 

^ = Vx^xfi). (27) 

In the limit where V{ = identically, the magnetic field is constant in time, and the neutrals simply 
move across the field lines, which are fixed in space. In this case, the spatial gradient of the field 
is given by the condition of hydrostatic equilibrium, equation (|2~5|), and the induction equation 
is automatically satisfied. This simple situation still applies for finite Vi as long as Vi <C%. To 
quantify matters, we rewrite the induction equation heuristically as 



where is the magnetic flux, defined by 



— = 2irBr. (29) 
or 



Equation (25) can be rewritten as 



° _ 2 K^v n Br. (30) 
at V n 



- 9 - 



Let L be a characteristic length of the problem, and t n a characteristic time for the motion of the 
neutrals; then 

v n ~ — , (31) 
and to within order of magnitude equation ( pOj ) reads 

d$ Vi BL 2 v,- $ 

(32) 



dt v n t n v n t n 

Equation ( |32|) shows that the characteristic time for the magnetic flux to change is a factor 
v n /vi longer than the characteristic time for the motion of the neutrals. Therefore, as long as 
v% /v n <C 1 the induction equation need not be solved separately, and the strength of the magnetic 
field is given by the condition of hydrostatic equilibrium. This means that the magnetic and 
hydrodynamic parts of the problem decouple, resulting in an enormous simplification. 

The validity of the condition Vi/v n <C 1 can be checked a posteriori by solving for Vi from 

d<5>/dt 



equation ( |28|) 



once B and are found, respectively, from eqs. (|25|) and (pS|). 



(33) 



Finally, the assumption that the thermal stresses are negligible compared to the magnetic 
stresses can be checked a posteriori by evaluating a modified form of the parameter [3 in 
equation (g): 



2.2. Dimensionless Expressions 

It is convenient to cast all the equations in dimensionless form. Let 

p n (t,M) = p C)0 p(r,m), (35) 

r(t,M) = R f(T,m), (36) 

and 

B (t, M) = B b (r, m) , (37) 

where p c> o, Rq, and Bq are, respectively, the initial central density, the initial radius of the cloud, 
and the magnetic field at the edge of the cloud, and where we have introduced a dimensionless 
time r and a dimensionless mass m: 



t = t AD T, 

M = Mq m. 



(38) 
(39) 
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Here Mq is the total mass of the cloud, and £ad is a characteristic ambipolar diffusion timescale, 

tAD-^^ = 4.96xl0 6 (-^) _1/2 yr, (40) 
47rG/9 V 2 V10 3 cm 6 J 

where n C) o is the initial number density at the center of the cloud. Note that the ratio tAo/ts, 
where tg is the characteristic free-fall time for the inital state, is given by 

= 147-^- =4 7 (41) 

or, in terms of the quantity vg used by Mouschovias and collaborators, 



t s V37r 2 ; 



is s . (42) 



Let us define the parameter e, which is inversely related to the initial degree of concentration 
of the cloud: 



-wf9 L - = ^- (43) 



_ M (p) 

-Ro 3 Pc,0 Pcfi 

Also, let Vi and v n be, respectively, the dimensionless ion and neutral velocities, defined to be 
positive, 

V n = -V n ,0 V n (44) 

and 

Vi = -VnflVi. (45) 

Here, 

_ GM R 

Vn <° = u n 2 1/2 = e T+ — ( 46 ) 

KniRo Pc,0 1 3tAD 

is a characteristic velocity of the neutrals. Note that the initial velocity of the neutrals at the edge 
of the cloud is v n = v n ,o[pc,o/ Po{Ro)] 1/2 > 1- 

To characterize the magnetic field, we define a parameter S as the initial ratio of gravitational 
to magnetic pressure at the edge of the cloud: 

The magnetic field can also be characterized in terms of the magnetic critical mass M$ , defined by 

M * = ^2 $ ' ( 48 ) 

where c$ is a numerical constant that depends on the distribution of the mass to flux ratio in 
the cloud; for a uniform spherical cloud with a constant magnetic field, c$ ~ 0.12 (Tomisaka et 
al. 1988b). Clouds with M > M$ (magnetically supercritical) are dynamically unstable and will 
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collapse even if the magnetic field is perfectly coupled to the neutral material. In terms of the 
characteristic value for the initial magnetic critical mass M<j>o = c^ttR^Bq/G 1 / 2 , we have 

M$Q TTC$y/2 



M 



S 1 / 2 



0.53 

W 2 ' 



(49) 



We shall focus on clouds that have small values of 8, corresponding to magnetically subcritical 
clouds {M§q/Mq > 1) in which magnetic forces are initially stronger than gravitational forces. 



Using the definitions in eqs.(|35|)-(f4"7|), our basic equations 



Dp 
Dt 

df 3 
dm 

dtf_ 

dm 







dm \ n 1 / 2 



£ 
P 



r m 



and equation ( |2l|) now reads 



m 



r 2 p x ' 2 ' 



), and (25) take the form 
(50) 

(51) 
(52) 

(53) 



The solution of equation ( p0| ) requires the specification of an initial condition p = p (r = 0, m) 
In addition, equations (^TJ) and (|52[) are subject to the boundary conditions 



f (r 



0, m 
df 
d7 



1) 



m=0 



1: 
0. 



and 



b (r, m 

respectively. The boundary conditions in eqs. 



D = h 



(54) 
(55) 

(56) 

are obvious, while the boundary condition 
given by equation (56) means that the magnetic field at the cloud's edge is constant. This 
boundary condition is appropriate for a dense cloud surrounded by a low density medium threaded 
by a field Bq. 

As an a posteriori check, the assumption of quasi-static motion can be quantified with the 
help of the quantity a defined in equation (]26f ) 

4ttG f 2 D v„ 



a 



k n 2 m Dt 
1.2 x 10~ 2 



f 2 Dv r 



m Dt 

where we used the value of k n i in equation (JTl|); the ratio Vi/v n is given by 

Vi 3 1 d(p/dT 



€ V r . 



(57) 



(58) 
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where <fi is the dimensionless magnetic flux defined by 

$ = ttB Rl <P; (59) 

and, finally, by defining a central value of the parameter (3 in equation (||), we can quantify the 
role of thermal stresses with 



0' 



87t/cb n C) o T dp /dm 
Bq db 2 /dm 
dp/dm 

Po p. u2 IPl , (60) 
db 1 [dm 



where n Ci o is the initial number density at the center of the cloud. 

To summarize, by putting our problem in dimensionless form, we have reduced the dependence 
of the solution to only two free parameters, e and 5, together with a specification of the initial 
mass distribution in the cloud. The parameter e is the inverse of the initial degree of concentration 
of the cloud (eq. fjSf ), and the evolution of the hydrodynamical variables f, v, and p depends 



only on this parameter. The second free parameter, 5 (eq. [47]), is the initial ratio of gravitational 
to magnetic pressure at the edge of the cloud, and parameterizes the initial mass-to-flux ratio of 
the cloud (eq. ^9|). The evolution of the magnetic field and the ion velocity depends on 5 and e. 
Also, we will show below (eq. |83|| ) that the value of 5 determines the regime where our static-ion 
approximation is valid. To convert our dimensionless solution into physical units, it is necessary to 
specify the cloud mass Mq, the initial cloud radius Ro, the magnetic field at the edge of the cloud 
Bq, and the ambipolar diffusion time £ad (eq. 0]) or the initial density at the center of the cloud 

Pc,0- 



2.3. A General Solution 



The key to the solution of the system of equations (pQ)- (|53|) is the solution of equation (|5Q). 
This can easily be achieved by defining the new dependent variable 

X = -^H (61) 

and the new independent variable 

f = i (62) 
m 

Using these transformations, equation (E0) takes the standard form 



a(x) § + ^ = °' (63) 

where 

o(x) = ~- (64) 
2x 
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Equation ( p3|) is a one-dimensional conservation law: if Xo{Q = x( r = 0)0 is the initial 
condition for \i then the solution of equation (]63|) is given by 



x(r,0=Xo[^-a(x)r]. (65) 



Two things are remarkable about the solution given by equation (65). First, it represents 
an implicit equation for x, and one hopes that for certain forms of the initial condition xo (0 
this equation can be solved algebraically. Second, note that the solution for \ has the form of a 
wave with velocity a(x), and therefore there exists the possibility of wave steepening and shock 
formation, i.e., the formation of a singularity. 



With a(x) given by equation (p4|), the solution in equation (|6q j becomes 

x(r,0 = Xom, (66) 

where 

V> = 1 + T -p 1/2 . (67) 
Because x = (p 1 ^ 2 ?) ; we have the solution for the density in terms of its initial distribution: 

p(r,e) = v 2 /5o m- 



To proceed, we must select a specific form for the inital density distribution /S (£)- A form 
that roughly approximates models for molecular cores is 

m p / i \ -p 

Po = 7- r-xp = 1 + — r > (69) 



(m c + m) p \ m c ^ 

with m c a constant. We shall concentrate on two cases: p = 2, which approximates the density 
profile of a non-singular isothermal cloud, and p = 0, which provides useful insights into the 
evolution of the central region. First, however, we consider some features of the general case. 



3. The Case of Arbitrary p 



Because 



rr 

m = Mq 1 / 47rr 2 pair, 
Jo 



(70) 



equation ( p9| ) is an integral equation for p (f), the initial density as a function of radius. This 
equation is readily solved to give 



Po 



1 (p + 1) r 



JT 

p+1 



em. 



(71) 
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Inserting this result back into equation ( p9[ ) gives 

(p + 1) r 



m 



1 



3n4r 

p+i 



em r 



1 



(r = 0), 



(72) 



and because, initially, m = 1 at f = 1, this result provides the relation between e and m c for our 
model, 

p + 1 



1 y +1 

1 + — I -I 



(73) 



We now evaluate the solution given by equation (6£) for the initial density distribution we 
have adopted, 

p = iP 2 (l + —Y P . (74) 



m c ip ) 

Using equation ((6?]) for we obtain the general solution 



m 
m r 



1 + Jp 1/2 



1 T 

+ 



2/p 



5V2 2 



(75) 



We obtain the central density p c by setting m = 0, 

4 



(76) 



(2-r) 2 ' 

Thus, the solution becomes singular at r = 2. Near the center of the cloud, the density becomes 



very large as r — > 2. We can expand equation (75) as this limit is approached and find 



m 



1 + 



2 + p 



2p/5 



1/2 



(77) 



We conclude that at r = 2 a core of mass (2/p)m c collapses to a singularity. This behavior is 
an artifact of our neglect of thermal pressure. Had we included thermal pressure, the singularity 
would have been avoided, and instead a high density object — a protostar — would have formed at 
r ss 2. 

One other feature of the solution at r = 2 is worth noting. Equation ([77]) shows that the 
mass just outside the singularity scales as /P 1 / 2 . Assuming a power law for for /5(f) there, we find 
p(r = 2, f ) oc f~ 2 : at the instant the singularity forms, the density just outside has the same radial 
dependence as a singular isothermal sphere. However, this gas is not static. 



The conditions for the solution found here to be valid are more easily discussed for specific 
models, so we now consider the collapse of a homogeneous cloud. 
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4. The Homogeneous Cloud 

The case of an initially homogeneous cloud (p = 0, e = 1) is of particular interest because a 
complete solution can be obtained analytically, and because it gives some insights into how the 
collapse of the inner core of a stratified cloud takes place. Here we have po = 1 an d f> = tp 2 , so 
that the complete solution is given by 



and 



f 3 



2 ' 
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(78) 
(79) 
(80) 

(81) 



1 3/2 




(82) 



There are several things that are worth noting. First, p remains uniform even at late 
times, just as in the collapse of a non-magnetic spherical cloud (Spitzer, 1978). The collapse is 
homologous, with v oc r, since the acceleration rises linearly with r for a homogeneous cloud. 
Figure 1 is a plot of r as a function of r for various values of m, in a cloud with 5 = 4 x 1CT 3 , 
which is appropriate for a cloud of Mo = 5M , Rq = 0.45pc, and Bo = 30 pG; such a cloud has 
M$o/Mo = 8.3, which is quite magnetically subcritical. Next, note that the expressions for the 
magnetic field and the flux introduce a characteristic time, r cr , which we choose to define as 

2-r cr = 2(3<5) 3 / 8 . (83) 

For the case at hand, r cr = 1.62. 

Because the outer parts of the cloud accelerate first, one is tempted to conclude that it is 
the cloud as a whole that first becomes gravitationally unstable, rather than the central region. 
This hypothesis can be tested by evaluating M$/M as a function of time and mass. For flattened 
clouds, it is customary to evaluate this ratio for a flux tube. In our spherical geometry, however, 
we evaluate M(r), the mass inside a sphere of radius r, and compare this with the flux 3>(r). Had 
we included the pressure that makes the cloud spherical, the condition for gravitational instability 
would have been that M(r) exceed both the Jeans mass and the magnetic critical mass M$. 
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Having M(r) > M$(r) is thus a necessary condition for gravitational instability. Figures 2a and 
2b show 

M& M$ /0\ 0.53 



m m w/ 51/2 



rn 



(84) 



where we have made use of equation (f49[). Note that at r = the cloud is magnetically subcritical 
everywhere (M$/M > 1), and that the outer edge of the cloud first becomes unstable at 
t ~ 1.62 = r cr . More precisely, at r = r cr , we have 4>{m = 1) = 1.935 1 ' 2 , so that M^/Mq = 1.03, 
independent of the value of 5. Thus, r cr is very nearly equal to the time at which the cloud first 
becomes supercritical. 

These results show that in a homogeneous, cold, and magnetized cloud the collapse due to 
ambipolar diffusion proceeds from outside-in 0. The inner regions of the cloud are gravitationally 
stable, and collapse only because of the weight of the overlying material. By contrast, the collapse 
of a non- magnetized singular isothermal sphere occurs inside-out (Shu 1977): such a sphere is 
unstable at all points, but the collapse begins at the origin, where the free-fall time is the shortest. 

Finally, we determine (frcore, the residual flux trapped in the cloud when it collapses to 
a singularity at r = 2. Taking the limit of <fi(m = 1) as r — ► 2 in equation (|82[), we find 
0core = (2<5/3) 1//2 . As a result, we have 

-MT = {-W)^ e = 0A3 - (85) 

Our model thus shows that the trapped flux is sufficiently small that the collapsed core is 
magnetically supercritical, as expected. 

Are these results consistent with our assumptions? The validity of the quasi-static 
approximation can be readily assessed with the help of the parameter a — the acceleration of the 
neutrals in units of the local gravitational acceleration (eq. p6[ ) — and the expressions in eqs. 
and fl80|), which yields 

r 2 Dv n 1 

(86) 



m Dt 6 
Therefore 

a = 2.1 x 10" 3 , (87) 

and the approximation of quasi-static evolution is very good, at least so long as our static ion 
approximation holds. 

Next, the assumption of negligible ion velocity compared to the neutral velocity has to be 
checked. Figure 3 is a plot of Vi/v n as a function of the Lagrangian coordinate M (in units of 



4 It is interesting to note that this outside-in collapse has been found by Goldreich & Weber (1980) in a completely 
different context. They studied the core collapse in supernova progenitors with an equation of state P oc p 4 ^ 3 . In 
our case, this equation of state becomes appropriate once the magnetic field has been amplified considerably, so 
that b S> 1 (see eqs. j78| and |^l[). Goldreich & Weber's results show that the dynamical evolution of the core is 
homologous, with the outside of the core accelerating first. 
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the cloud's mass Mo) for four different times. It can be seen that up to t/t ad = 1-62 — i.e., at 
r cr , when the cloud as a whole becomes magnetically supercritical — Vi/v n <^ 0.3; therefore our 
assumption of negligible is justified a posteriori, up to the point where the cloud becomes fully 
supercritical. Once the cloud is supercritical, V{ begins to approach v n , and our solution breaks 
down. This accounts for the late-time behavior of M$/Mo seen in Figure 2b: after r = 1.62 this 
ratio does not change much, as expected when the magnetic flux starts to become frozen onto 
the neutrals. We conclude that the static ion approximation breaks down at about the same 
time the cloud becomes magnetically supercritical; this ultimately leads to the breakdown of the 
quasi-static approximation for the neutrals. 

We can analytically verify that the static ion approximation breaks down when the cloud 
becomes supercritical. Note that the maximum value of Vi/v n occurs at the center of the cloud. 
According to equation (^), Vi/v n depends on the ratio of dcft/dr to dcj)/dr. Near the center of the 
cloud, we have 4> ~ irf 2 b (to <C m c ). Evaluating the partial derivatives, we find 

(<) = i + 2 fl^r 

\ v nJ m=0 V 2 — T cr / 

This result confirms that the ions are approximately static prior to r cr , when Vi/v n = 1/3. Note 
that this estimate indicates that the ions and neutrals become comoving as r — > 2, as they do in 
reality. 

A remarkable feature of the solution discussed in this section is that it is completely 
independent of one of the basic parameters of the problem, the ratio of the ambipolar diffusion 
time to the free fall time {vg oc £ad/%)- The reason for this is now clear: our solution is valid 
only during the magnetically subcritical, quasi-static stage of evolution, when the timescale is 
determined by £ad alone. The solution sets the initial conditions for the subsequent supercritical, 
dynamical stage of evolution, when the timescale is set by tg. 



5. A Stratified, Non-Singular Cloud 



We now consider a more realistic case, namely, that of a stratified, non-singular cloud. For an 
initial density profile of the form 

2 

p(r = 0,m) = - ^— j (89) 

(m c + to) 

(i.e., choosing p = 2 in eq. H) equation © can be inverted algebraically, and we obtain 

4 
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1 + 



T + \ 1 + 



2t- 



2 ' 



(90) 



Recall that this case is a good approximation to a stable, non-singular, isothermal sphere, which 
is relevant to the study of real clouds. 
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Armed with the solution in equation (|90|), we can solve for r (r, m) using equation ( |5l| ) and 
the boundary conditions^] in equations (E3)-(|53): 



r ( r, m) 
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(91) 



What is the nature of the solution given by equations ( |9C|) and (91)? Recall that at the time 
r = 2, a core of mass M c = m c Mq collapses into the origin, leaving behind an envelope with mass 
(1 — m c )Mo. Evaluating p and f 3 at r = 2, they take the form 



p(r = 2,m) =4 
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(92) 



(93) 



respectively. Note that for m < m c the expression between the square brackets in equation (p3|) 
vanishes, and f = identically. We conclude that in a cloud that is initially everywhere 
magnetically subcritical, a core with a fraction m c of the total cloud mass will, through ambipolar 
diffusion, eventually become supercritical, and collapse into the origin. The relationship between 
m c and the initial degree of concentration of the cloud, e, given by equation (|73|) with p = 2, is 
plotted in Figure 4. We find that m c = 1 for e = 3/7, and therefore clouds with e > 3/7 will 
collapse as a whole, leaving no envelope behind. 



5.1. Evolution Prior to Collapse 

Further insight into the nature of the solution can be gained by examining in detail the 
time evolution. Here we consider the evolution up to the formation of a singularity in the 
origin of a cloud with e = 0.1 (corresponding to m c = 0.26) and 5 = 4 x 10 -3 (or, equivalently, 
M^q/Mq = 8.4). These values correspond to a cloud of mass Mq = 5M0, initial radius 
Rq = 0.45pc, a magnetic field strength at the edge of the cloud B$ = 30 pG, and an initial particle 
density at the center of n C) o = 2.2 x 10 3 cm -3 — yielding a characteristic ambipolar diffusion time 
t AD = 3.30 Myr. 



5 The boundary condition given by equation (£^) is strictly valid for r < 2 because (r = 2, m < m c ) is a singular 
point of equation ([Hi]). However, the solution for f(r, m) given by equation can still be used for r > 2 and 
m > m c . The reason is that the evolution of shells enclosing m > m c depends, by Newton's theorem, on m only and 
not its distribution — whether pointlike or distributed over a finite volume. Therefore, the neutral's velocity v n and 
the location of each shell f are smooth functions of r across r = 2 for m > m c . 
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Figure 5 is a plot of the Lagrangian position r, neutral velocity v n , density p, and magnetic 
field strength B. Focusing on Figure 5a, it is noteworthy that by r = 1.86 the radius for m = 1 
is reduced by a factor ~ 0.6 from its initial value , while for m <^ m c this change is by a factor of 
order ~ 0.2. Moreover, note how the change in radius is roughly the same for all m in the region 
m <^ m c , i.e, the evolution of the core resembles that of a homogeneous cloud. This result is 
further illuminated by considering the behavior of v n . Initially the absolute value of v n is largest 
at the outer edge of the cloud, but by t/t ad ~ 1.6 a maximum (in absolute value) starts to develop 
at m ~ m c , which becomes more accentuated as time goes on. In other words, the dynamical 
behavior of the core with m <J m c resembles that of a homogeneous cloud, and its collapse also 
proceeds from outside-in. 

Finally, it is interesting to examine Figure 5b and note that by t/t ad = 1-8 the central density 
has increased by a factor ~ 100, while B at the center has increased by only a factor <^ 4 — as 
expected in a cloud whose evolution is dominated by ambipolar diffusion. 



5.2. Self-Consistency 

In this section we check the self-consistency of our results for our fiducial stratified cloud, 
namely, the validity of the quasi-static approximation, the neglect of thermal stresses, and the 
assumption that Vi <C v n , which allows us to neglect the induction equation. 



5.2.1. Quasi-static Evolution 

Figure 6 is a graph of a, the acceleration of the neutrals in units of the local gravitational 
acceleration, as a function of m and r. Note that, up to the time shown (r = 1.996), a <C 1, i.e., 
the acceleration of the neutrals is a very small fraction of the local gravitational acceleration, and 
therefore the assumption of quasi-static evolution is very good. Also, note that for m ^ 0.2, a is 
roughly constant, as in a homogeneous cloud. However, just as in the case of the homogeneous 
cloud, it is the breakdown of the static ion approximation that occurs first, and that ultimately 
leads to dynamical evolution. 



5.2.2. Static-Ion Approximation 

A key assumption in our approach is that the ion velocity is much smaller than the neutral 
velocity, so that the hydrodynamic and magnetic problems are decoupled. This key assumption 
has to be verified a posteriori by computing the magnetic flux $ (eq. f29|| ) and computing v% 



as prescribed by equation (33). The evolution of the magnetic flux as a function of mass and 
of position is shown in Fugure 7. Figure 8 is a plot of Vi/v n as a function of the dimensionless 
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Lagrangian coordinate M/Mq for five different times (in units of the characteristic ambipolar 
diffusion time £ad)- Note that for r <; 1, Vi/v n <C 1 obtains everywhere in the cloud, justifying 
our initial assumption. By r = 1.64, vi/v n ~ 0.35 in the inner few percent (by mass) of the cloud. 
This result is similar to the one we obtained for the homogeneous cloud above. By plotting M$/M 
(Figure 9), we see that a part of the cloud first becomes magnetically supercritical shortly before 
r = 1.64. Hence, just as in the homogeneous case, the static ion approximation breaks down when 
the cloud becomes supercritical (although in this case only the core region is supercritical). 

5. 2. 3. The Validity of the T = Approximation 

One of our main simplifying assumptions is that thermal stresses are negligible compared to 
magnetic stresses. The validity of this approximation can be checked a posteriori by looking at 
the value of /?', i.e, the ratio of the thermal to magnetic pressure gradients (eq. [Q|). Figure 10 
is a plot of (3' as a function of m and r for our fiducial cloud. Note that /?' <C 1 throughout the 
evolution; therefore, our assumption that thermal stresses can be neglected is justified. 

At a given time (e.g., for i/tAD = 0), our model predicts that the importance of thermal 
support would increase as one moves towards the center of the cloud until m <^ m c , where it starts 
to decrease. In a real cloud, one expects that, at the center, thermal pressure becomes important 
for a mass of order the Jean mass at the appropriate central density. Thus, one expects to 
increase towards the center of the cloud; whether at the center (5' ~ 1 depends on the particular 
magnetic configuration. Also, as ambipolar diffusion proceeds, the central density increases, and 
the mass — and size — of the region where thermal stresses are important decreases with time (recall 
that the Jeans mass oc (n c ) _1//2 , where n c is the central density). In our case, the local maximum 
of (3' obtains at a fixed value of m — namely, m = m c — because of our choice of initial conditions 
and our neglect of thermal stresses. However, this is not critical as long as we find that j3' <C 1 
throughout the cloud, because, as mentioned above, the Jeans mass decreases with increasing 
density. Thus, at worst, our calculation is in error in a region that becomes increasingly smaller 
with time. 

As time increases, thermal pressure forces at the center may begin to dominate the magnetic 
stress due to the loss of magnetic flux. The answer to this question, again, depends on the 
particular magnetic configuration. In our case, we find that thermal stresses never become 
important during the quasistatic phase; this finding is in agreement with the results of Fiedler 
& Mouschovias (1992), who find that even during the dynamical collapse of the core, magnetic 
stresses dominate over thermal pressure at least up to a density of ~ 10 9 cm~ 3 . In our model, if 
there is no core to begin with-the constant density case- one never forms; (on the other hand, cores 
can form during the subsonic evolution of clouds supported by gas pressure — see Bodenheimer & 
Sweigert, 1968). Once a core has formed, then we find that the subsequent contraction of the 
cloud results in a decrease in the relative importance of thermal stresses, as shown in Figure 10. 
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Note that this predominance of magnetic support over thermal support does not contradict 
the observational evidence that the linewidths of dense cores are essentially thermal (see, e.g., 
Caselli & Myers, 1995, and references therein). Recall that, in our picture, the support of the 
cloud is due to quasistatic magnetic fields — as opposed to magnetic turbulence — and that the drift 
velocities during the quasistatic phase are an order of magnitude below the speed of sound (see 
Figure 5a). Therefore, in our model, the main contribution to the line widths is thermal motion. 



5.3. Comparison with the Numerical Work of Fiedler & Mouschovias 

Are our results relevant to the evolution of real clouds? To date, one can obtain an answer to 
this question only by comparison with state-of-the-art numerical calculations. Here we compare 
our results against those of Fiedler & Mouschovias (1992) for their baseline model, Model 1. To 
effect the comparison, we need to specify five parameters: the time to i n Fiedler & Mouschovias's 
results that we identify with r = 0, the corresponding central density p Cj o, the mass and initial 
radius of the cloud, Mq and Ro respectively, and the magnetic field at the edge of the cloud, 
Bq- The values of these last three parameters are fixed by Fiedler & Mouschovias's choices, 
namely Mq = 45.5 M , Rq = 0.75 pc, and Bq = 30 pG; the values of to and p C) o require further 
consideration. 

Because Fiedler & Mouschovias start their calculations with a homogeneous, uniformly 
magnetized cloud that relaxes to a quasi-static state, the only practical way to determine to 
and p c fl is by fitting their results for the central density during the quasi-static phase with our 
expression for the central density (eq. |7(| in dimensional form): 



p c (t) = 4 p cfl 



^Gp l Jht-t ) 



k,, 



(94) 



Using the values adopted by Fielder & Mouschovias in their Model 1 for the mean molecular weight, 



Ki, and {av} ni , we find that there is a continuum of equally good fits with n Cj o ^ 1.5 x 10 3 cm -3 — 
or, equivalently, e <^ 0.3 — and to ^ 8.0 Myr. However, we show below that geometrical constraints 
narrow the range to e ^ 0.1 (n Cj o <; 4.5 x 10 3 cm -3 ) and to ^ 10.0 Myr. We decide to adopt our 
fiducial value of e = 0.1, and a very good fit obtains with to = 11.85 Myr. Given that Mo and 
Rq are fixed, this choice of e translates into an initial central density n C) o = 4.4 x 10 3 cm -3 , and 
the characteristic ambipolar diffusion time for the cloud is £ad = 2.33 Myr. Finally, the choice 
-Bo = 30 pG translates into 5 = 4.1 x 10 -2 — a factor of ten larger than our fiducial value — with 
M^o/Mo = 2.6. 

Figure 11 is a comparison of the evolution of the central density and central magnetic 
field — this last obtained from evaluating B at m = 10 -3 — between our model and Fiedler & 
Mouschovias's. Focusing first on the central density, we see that even though to = 11.85 Myr, 
the fit is good starting at t ~ 8Myr and remains good until t ~ 16Myr (t/t ad ~ 1-8). After 
t = 16 Myr the two models diverge, and no fit can be found for any choice of to and n c o- This 
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is due to the onset of dynamical collapse in the Fiedler & Mouschovias calculation, while our 
model still evolves quasi-statically up to t ~ 16.5 Myr (refer back to Fig. 6). The reason for this 
discrepancy can be found in the evolution of B c in Figure 11: note how the strength of the central 
magnetic field increases faster in our model than in Fiedler & Mouschovias' — by t ~ 15 Myr the 

2 /3 

difference is a factor ~ 2, and thereafter we find that B c oc n c , i.e. flux-freezing obtains in 
the core. This larger magnetic field provides enhanced support against gravity and the onset of 
dynamical collapse is retarded. However it is interesting to note that at n c ~ 10 5 cm -3 , where our 
solution starts to differ significantly from that of Fiedler & Mouschovias, the free-fall time in the 
core is ~ 0.1 Myr and the ambipolar diffusion time is ~ 0.5 Myr. Therefore, our solution for n c 
agrees with that of Fiedler & Mouschovias up to the point where the evolution of the core becomes 
dynamical and the core is within a few free-fall times from reaching the origin, notwithstanding 
the large discrepancy in B c . This is so even though the static ion approximation breaks down 
earlier than in our fiducial cloud — we find that in the core Vi/v n ~ 1 by r ~ 1.4 due to the fact 
that 5 is an order of magnitude larger than the fiducial value used in our numerical examples. 

So far we have only proven that our solution compares well with that of Fiedler & Mouschovias 
in the central region. Given the discrepancy in the late-time evolution, it is interesting to compare 
the evolution in other parts of the cloud. Because of the differences between the two calculations, 
the only meaningful comparison is that between our results and those of Fiedler & Mouschovias 
for the midplane density of their cloud; this comparison is presented in Figure 12, which is a 
plot of density as a function of radius for three different times, t = 10.2 Myr, t = 15.1 Myr, and 
t = 16.0 Myr. Two things are noteworthy in this figure. First, even though the central density 
disagrees by a factor of ~ 2, the overall profile is well reproduced^. Second, note how the density 
profile in the envelope is also well reproduced (p oc r -2 ), but Fiedler & Mouschovias's values 
are consistently larger than ours (again, by ~ 2); this is a natural consequence of the different 
geometries. Although Fiedler & Mouschovias's model is fairly spherical — with an aspect ratio 1:3 
in the envelope — our cloud is strictly spherical, so our density in the outer parts has to be lower 
for the two clouds to have the same total mass. This constraint guided us in choosing e = 0.1 for 
the comparison. Otherwise, for e > 0.1, the spatial extent of our core and our envelope density 
would be larger than those in Fiedler & Mouschovias's cloud, and the clouds' masses would not 
match. 

These results confirm that the assumptions of negligible thermal stresses and of static ions 
have minimal consequences for our results for the hydrodynamical variables until the latest stages 
of gravitational collapse. 



However, note that for t = 16.0 Myr, where the central density is reproduced most accurately, the drop in density 
at the edge of the core is steeper in the Fiedler & Mouschovias model than in ours. This is a natural result of 
neglecting thermal pressure: the density contrast between the center and the edge of a stable, thermally supported 
cloud is 14 : 1 (Bonnor, 1956), and the mass of the core decreases with increasing density as p -1 ^ 2 ; in our model the 
mass of the core is fixed, while the density contrast between the center and the edge of the core increases with time 
from 4 : 1 at r = (eq. @). 
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5.4. 



The Initial Protostellar Mass 



What fraction of the cloud becomes the actual protostar? In our model, the answer is m c Mo, 
which is fixed ab initio once the initial degree of concentration of the cloud e is fixed. This 
characteristic mass corresponds to the inner fraction of the cloud that first becomes magnetically 
supercritical and collapses, detaching itself from the remaining envelope. However, m c is a free 
parameter in our model; what determines in nature, then, the protostellar mass? Notwithstanding 
the shortcomings of our model, here we show that it provides insight to this question. 

As we noted before, the stability of a cold cloud is governed by the ratio M$/M, where M$ 
is the magnetic critical mass (eq. pjfl ). Those regions of the cloud with M > M<j> are dynamically 
unstable and will collapse even if the magnetic field is perfectly coupled to the neutral material. 
Figure 9 is a plot of M$ /M as a function of time and the Lagrangian coordinate M for our fiducial 
cloud (e = 0.1 and 5 = 4 x 10~ 3 , with M^^/Mq = 8.4). Two things are noteworthy in this figure: 
first, the cloud is initially magnetically subcritical throughout, and as time increases the region 
M ~ M c first becomes supercritical at t/iAD ~ 1-5. Second, and more important, at any given 
time the ratio M$/M has a global minimum as a function of M, and this minimum occurs at 
finite M. Thus, a region of finite mass first becomes magnetically supercritical, and the collapse 
proceeds inward from this region. This result is readily understood: for small enough values of 
m, both B and p are fairly constant, so that M$ oc BR 2 oc Af 2 / 3 . On the other hand, in the 
outer parts of the cloud p oc r~ 2 and B ~constant, so that Af$ oc R 2 oc Af 2 . Thus altogether 
M$/M oc M -1 / 3 + KM, where K is a constant, and this expression has a minimum at a mass at 
which the two terms are comparable. 

Is this result relevant to clouds with a finite temperature? The answer is yes — as long as 
magnetic stresses dominate outside the thermally supported core. In a real cloud, there are two 
characteristic masses, namely a Jeans mass which can be supported by thermal stresses only, 
Afj, and the magnetic critical mass M$; and both are time dependent. Therefore, at any time 
the stability of the cloud is now determined by the ratio (Af$ + Mj)/M (McKee, 1989). If the 
core is thermally supported, this function will also have a global minimum as a function of Af. 
Just as before, this minimum will occur at a mass at which the two terms are comparable — i.e., 
M<j, ~ Mj — so that Af is of order Afj. We conclude that the characteristic protostellar mass is of 
order the Jean mass at the onset of dynamical evolution (see also Larson 1985, 1995): 



where c s is the speed of sound and Pq is the thermal pressure confining the cloud, which is typically 
somewhat greater than the thermal pressure in the diffuse ISM (Bertoldi & McKee 1992). Note 
that if the core is thermally supported, then the thermal pressure is comparable to the magnetic 
pressure and the characteristic stellar mass Mj oc 1/B (cf. Basu & Mouschovias 1995b). In our 
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(95) 
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model, the role of the Jeans mass is taken by the core mass M c . 



5.5. The Protostellar Magnetic Flux 



Our model allows us to estimate the initial magnetic flux trapped in the protostellar core, 
&core(T = 2). We have already found the trapped flux analytically for the homogeneous case; 
here we determine it for the p = 2 case. We can do this despite the breakdown of the static ion 
approximation before r = 2 because the breakdown of the approximation leads to V{ ~ v n , so 
that flux freezing obtains. As a result, the value of § COTe (T = 2) is essentially the value of $ cor e 
at the time the approximation breaks down. This is shown in Figure 13, which is a plot of the 
evolution of the magnetic flux in the core, ^coro in units of the cloud's initial magnetic flux, <&q, 
for S = 10 -3 and different values of e. Note that the amount of flux trapped in the core at r = 2 
is a decreasing function of e, or, equivalently, of m c . 

How does < ^corc('7" = 2)/3>o scale with m c and <5? We can answer this question by examining 



equation (52). To isolate the dependence on m c , let us introduce a new Lagrangian variable 
y = m/m c . Then, equation (52) reads 



dy 



-5m c 2 ^j. 
r 



(96) 



According to equation (|9l|), we can write r 4 = (m c e) 4 / 3 f(r, y). Let us assume m c <C 1, so that 
e ~ 3m 2 (eq. |F3] with p = 2). We then obtain from equation (|96|) 



6 2 (r,y) = l + 



35 [ x l m * y dy 



mf. 



f(r,yY 



(97) 



In the core (y < 1), and for r > 1, we have b 3> 1 (see Figure 5b), and therefore the first term on 
the r.h.s of equation (|97D can be neglected. Furthermore, because the gradient of b is largest just 
outside the core (again, see Figure 5b), the integral in equation ( |97| ) depends weakly on the upper 
limit of integration, so that 



6(t > l,y < 1) oc 5 1/2 



in 



-i 



(98) 



As for the magnetic flux in the core, by rewriting equation ( |29| ) in Lagrangian coordinates 
and integrating we obtain, in units of <I>o> 



'core — 

3 Jo rp 



1 bdy 



(99) 



Because p depends on r and y only (eq. |9^j) whereas f scales as (m^) 1 / 3 , we finally obtain that 
at late times 

</w(r > 1) ocm c 5 1/2 . (100) 
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The scaling in equation (100) describes our numerical results quite well, and we find the coefficient 
of proportionality to be 1.16. By analogy with equation (^), we find that the value of M$/M in 
the core is 

^! = W^ = 0,61, (101) 

independent of m c or 5. Just as in the case of the homogeneous cloud, the flux is reduced until the 
core is supercritical; the more subcritical the core is to begin with, the greater the flux lost. 



5.6. Post-Collapse Evolution 

A particularly attractive feature of our solution is that we can follow, with the help of 
equation (p0|), the evolution of the envelope with m > m c after the collapse of the core at r = 2. 



As in the pre-collapse evolution, one should check for self-consistency regarding our 
quasi-static approximation, the neglect of thermal stresses, and the assumption that vi <C v n . 
The validity of the quasi-static approximation is assessed with the help of a (eq. J26|), the local 
ratio of the neutral's acceleration to the gravitational acceleration; similarly, the validity of the 
T = approximation is verified by evaluating j3' (eq. |6(|). On the other hand, the assumption of 
negligible motion of the ions requires special consideration. 

Verification of the assumption Vi/v n <C 1 requires $(t,r) to be known. It is difficult to 



evaluate $(t, r) after the collapse of the core in the same way as for r < 2 (integrating eq. [29| 
outwards from r = 0) because, as more and more mass shells accrete onto the protostar, the 
protostellar magnetic flux increases with time, and its evaluation requires the solution of a partial 
differential equation. However, one can evaluate the magnetic-flux loss at the edge of the cloud and 
find <I>(Af = Mo) as a function of time in the following way. Adding and substracting v n d&/dr 
to the l.h.s of equation (p§|), substituting d<&/dr from equation (29), and switching to Lagrangian 
coordinates we obtain 

2TTv n [ 1 - — ) Br ~ 2irv n Br, (102) 



Dt V v u . 

where the last step follows for Vi <C v n . Note that at the edge of the cloud we have B = Bq, 
so that equation (102) may be integrated to evolve the magnetic flux at M = Mq. With this 



value at hand, we can then integrate equation (^9|) inwards. It turns out that integrating inwards 
equation (|29|) with the initial value found from equation ( |102| ) starts to breakdown at the point 
where Vi ~ v n (we obtain negative values for $; recall that v n is negative definite), but, as long as 
a <C 1, the assumption <C v n is very good throughout the quasi-static envelope^. Finally, we 
find that the assumption of negligible thermal stresses is very good throughout the post-collapse 



7 The same method may be applied to find $ during the pre-collapse stage. We find that as long as Vi v n the 
two methods (integrating outwards from r — and integrating inwards from M = Mo) agree very well. However, 
as soon as approaches v n in the core, the integration from M = Mo breaks down for m < m c , although it yields 
Vi <S v n throughout the envelope. 
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evolution and that our results are only limited by the breakdown of the quasi-static approximation. 
We arbitrarily deem the evolution to become dynamical where a > 0.1. 

Figure 14 is a plot of the neutral velocity and density in our fiducial cloud as a function of 
position r for six different times after the collapse of the core; for each curve, the smallest value of 
r is that where a = 0.1, while the largest corresponds to the edge of the cloud. For Mq = 5 M and 
Ro = 0.45pc — and recalling that e = 0.1 — we obtain n Cj o = 2200 cm -3 and v h q = 4.5 x 10~ 3 kms , 
therefore the typical densities in the envelope are ;> 750 cm~ 3 and the maximum velocities are 
~ —0.10 kms -1 — a factor <^ 2 smaller than the speed of sound for T = 10 K. 

Two things are noteworthy in Figure 14. First, after the collapse of the core at t/txD = 2.0 
the density profile starts to flatten as the velocity profile steepens. This is a natural result of the 
increased acceleration of the neutrals. Second, the point where the evolution becomes dynamical 
moves outward with time, and therefore the collapse of the envelope is from inside-out after the 
collapse of the core (see also Foster & Chevalier, 1993). For our fiducial solution this collapse front 
moves at a speed ~ 0.05 kms -1 (i?o/0.5pc) (n Cj o/10 3 cm -3 ) 1 / 2 — smaller by ~ 4 than a typical 
speed of sound in dark clouds. 



5.7. Mass Accretion Rates 

We can use our results to evaluate the mass flow rates in regions where the flow is quasi-static, 
and to estimate the accretion rate onto the protostar after it forms. In spherical symmetry the 
mass accretion rate at any point r in the cloud is given by 

M = ATrr 2 v n p = M AD mp 1/2 , (103) 

where we have used the non-dimensional expressions in eqs. (^), (^), (fH|), and (|53|). The 
characteristic accretion rate due to ambipolar diffusion is 

M 



M AD 

*AD 



= 41x10 -m s ^(^) 3/2 (^)- 3/ V/ 2 (104) 

where we have used equations ([£(]) and ([R]); recall that e is the initial ratio of mean to central 
density. Figure 15 is a plot of M (in units of the fiducial mass accretion rate Mad) as a function 
of radius in our fiducial cloud for both the pre- and post-collapse stages. 

Focusing first on the pre-collapse evolution, we note that for our fiducial cloud (e = 0.1) and 
reasonable values of the cloud's mass (Mo ~ a few M Q ) and initial radius (Ro <; lpc) the value 
of M for most of the cloud is ~ few x 10 _7 MQyr _1 . However, at any given time t/t ad < 2, 
M has a maximum as a function of r; this maximal value increases with time, and its location 
moves inward with M(r = 0) — > oo as r — > 2. This feature is a direct result of our assumed initial 
condition, which results in a constant core mass as a function of time. 
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As for the post-collapse accretion rates, M(r) is plotted in Figure 15 up to the point where the 
quasi-static approximation breaks down (a > 0.1). The results in Figure 15 show, again, that to a 
good degree of approximation the mass accretion rate in the quasi-static envelope (t/Rq ^ 0.1) is 
constant — with M ~ few x 10 M^yr -1 . 

How do these results depend on the free parameters, namely Mq, Rq, and e? Using 
equation (|104|), equation (|103|) can be rewritten as 



M = 4.! x lO-M 0yr -> (^L) V2 (^)- S/2 ^. (105, 

We find that in the range 0.01 < e < 0.2 the combination e -1 / 2 mp 1 ^ 2 ~ 0.6-1.6 (depending on 
the values of r/R$ and t/iAD where the expression is evaluated); this translates — holding Mq 
and Rq constant — into a change <J 3 in the mass accretion rates. Therefore, we conclude that 
for reasonable choices of the initial degree of concentration (0.01 < e < 0.2), the cloud's mass 
(Mq ~ a few solar masses), and initial cloud radius (Rq = 0.5-1.0pc) the mass accretion rate in 
the quasi-static envelope is ~ few x lO^M^yr -1 . 

Finally, we address the question of the mass accretion rate onto the protostar, M*. By 
continuity, one would expect M* to be comparable to the mass accretion rate in the quasi-static 
envelope. We can answer this question more quantitatively with the help of our expression for 



r(t,M) (see eq. |9l[]). Setting f = in equation (91), we can solve for the instantaneous value 
of m at the origin, i.e., the protostellar mass as a function of timef]; and by computing its time 
derivative, obtain M*(t). Since we cannot treat the dynamical stage of evolution, this procedure 
will become reasonably accurate only after about one free fall time after the collapse. However, 
although we cannot treat the rate of accretion accurately during this stage, we can obtain a good 
estimate for the integrated accretion-i.e., the total mass of the protostar that forms during the 
dynamical phase. Figure 16 is a plot of M* and M* (in units of, respectively, Mq and Mad) as a 
function of time after the formation of the protostar for various values of e. 

In our model, the formation of the initial protostar corresponds to an infinite accretion 
rate. For r-2 < 1, we find that M^/Mq — 1 oc (r — 2) 1 / 2 , and therefore M — and the accretion 
luminosity — diverge as (r — 2) -1 / 2 (though this is an integrable singularity, and the total energy 
radiated by the shock vanishes as (r — 2) 1 / 2 ). However, this is an artifact due to neglecting 
thermal stresses. Recall that in our model the protostellar mass is a discontinous function of 
time near r = 2, with M* = as r — > 2~ and M* = m c Mo as r — > 2 + ; the infinite accretion 
rate is a result of this discontinous behavior. When thermal pressure is included, it provides a 
retarding force whose effect makes the protostellar mass a continous function of time at r = 2. 



8 For r > 2 — and as a result of ignoring thermal pressure — f — is a singular point of equation However, 
we can still find M,(r) from equation ( |9lj ) by setting f — because a real protostar will have a finite radius. For a 
typical protostellar radius 7?, ~ 5Rq and an initial cloud radius Ro ~ 0.5 pc, the lhs. of equation ( |9lj ) evaluated at 
the protostellar boundary is of order 10 -20 — or zero, for all practical purposes. 
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The calculations of Foster & Chevalier (1993) illuminate this point. They followed the collapse 
of a non-singular cloud when only thermal pressure opposes gravity; therefore, their results are 
directly applicable to our phase of dynamical infall near r = 2. They find that M is finite at the 
instant the protostar is formed, and it rises very sharply immediately thereafter. After this initial 
rise, which lasts for about a tenth of a free-fall time, the accretion rate declines in very much the 
same way as in our model. 

For our fiducial choice of parameters, a tenth of a free-fall time at the density of dynamical- 
infall onset corresponds to (r — 2) ~ 2 x 10 -3 , or t 10 4 yr after the formation of the initial 
protostar. Assuming a radius for the hydrostatic protostellar core R* = 5Rq, we find an accretion 
luminosity 

L^SOWe-V* r^Vf 1/2 f-^V S/2 f^V 1 . (106) 

° c V2xio- 3 y \5M Q J vo.5p C y \5r q J v ; 

For our fiducial model (e = 0.1 and m c = 0.26, which corresponds, for Mq = 5M , to an initial 
protostar of 1.3 Mq), we find that the accretion luminosity ~ 10 4 yr after the formation of the 
initial protostar is 49 L@. For an initial protostellar mass M* = 0.5 M Q (obtained with e = 0.025, 
which gives m c = 0.1), we find, for the same time after the initial collapse, an accretion luminosity 
~ 15 Lq. This last value is only a factor ~ 2 higher than those of typical Class sources 
(Gregersen et al., 1997). 

As shown in Figure 16, the mass accretion rate remains very high for about a free fall time 
(~ 0.2 £ad) after the collapse, and, for the range in e explored, the protostellar mass increases by 
a factor ~ 2 during this phase. Thereafter, M* is comparable to the mass accretion rates in the 
quasi-static envelope, as expected. This strong dependence of M* shortly after the formation of 
the protostar has also been found by Foster & Chevalier (1993) and Hunter (1977) for the case of 
non-magnetic collapse — a result that is at strong variance with the predictions from the singular 
isothermal sphere model (Shu 1977). The later stages of accretion are more in accord with the 
constant accretion rate expected in the singular isothermal sphere model; particularly — when one 
takes into account that Mad °c e^ 1 / 2 (eq. [104]) — M* is essentially independent of e during the 



post-collapse, quasi-static phase. This similarity with the singular isothermal sphere model merits 
further attention, and we address it in the next section 



5.7.1. The Post-Collapse Mass Accretion Rate at Late Times 
and the Singular Isothermal Sphere Model 

The standard paradigm for the gravitational collapse of an interstellar cloud is the inside-out 
collapse of a singular isothermal sphere (Shu 1977). The cloud is assumed to be initially in an 
unstable equilibrium in which gravity is balanced by the pressure of an isothermal gas so that 
p = c 2 /27rGV 2 , where c is the isothermal sound speed. Collapse begins at the origin, where 
the dynamical time is the smallest, and an expansion wave moves outward at c. At a time t 
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after collapse begins, a mass 2c 3 t/G has been swept up by the expansion wave; of this, a mass 
0.975c 3 i/G is in a condensed core at the origin. 

In light of our results for the post-collapse, quasi-static accretion rate — namely, that at late 
times after the collapse of the inner core of a non-singular stratified cloud the mass accretion 
rate onto the protostar approaches a constant value — we wish to compare our predicted accretion 
rates at late times with that predicted by the singular isothermal sphere model. To do this, we 
note that at r = 2 the density distribution of the envelope is given by p oc r -1 / 2 (see Figure 14), 
the same density distribution as the singular isothermal cloud. At late times (t» 1), the mass 
shells that still are in quasistatic equilibrium satisfy m 3> m c ; in this limit, equation ( f72|) for 
m(r) reduces to m = f, and the initial density distribution (eq. |2q| ) after core collapse becomes 
p(r = 0) = (m c /m) 2 . Rather than measuring time in units of the central ambipolar diffusion time 
*AD) we use now the ambipolar diffusion time at the surface, t' AB = tAD[Pco//°(0, l)] 1 ^ 2 = tAD/ m c- 
The density distribution becomes 

p(r', to) 4 . .„ 

^77rTr = T T—r-2 o m/2 12 ™>0, 107) 

p{0, 1) [to — T + (m z — 2mT') ] -' 2 \ z 

where r' = i/^AD = m c r an d t' = corresponds to the formation of a singularity at the origin. 
The radius of the mass shell to is 

13 3 1 

f 3 = -m 3 - -m 2 r' + -tot' 2 + -(to 2 - 2mT') 3/2 (m > to c ). (108) 

This equation indicates that the shell reaches the origin at a time r' cr = |to. By this time, our 
approximations for the dynamics have broken down, but this nonetheless defines a characteristic 
time scale for the evolution. At this time, equation ( |1Q7| ) shows that the density has increased by 
a factor of about 5. 

If we approximate the acceleration of the gas by t> 2 /r, the parameter a (eq. p6|) becomes 



4vrG 

a 



k 



in 



r(0, to) 



r(r', to) 



p(Q,m) 



_p(r',m) 



(109) 



As discussed in §5.2.1, a is initially very small; it will begin to approach unity (so that our 
quasi-static approximation breaks down) only when f 3 becomes significantly smaller than its 



initial value. Equation (108) shows that this will occur only when r' is quite close to t' cx . It is 



therefore self-consistent to use to use t' cv as an estimate for the evolution time for the shell. 

The time scale for evolution t' cv enables us to define the characteristic velocity with which 
ambipolar diffusion advances into the cloud, 

KM 9 R 
v ad = -r-n— = T ~p ' (HO) 

r cr r AD ^ I AD 

which is constant. How does this compare with the expansion wave associated with the collapse 
of an unmagnetized, singular isothermal sphere, which moves outward at the sound speed c? 
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Defining an effective sound speed through the relation governing such a sphere, c 2 = 2irr 2 Gp, and 
using equation (40), we find 

3.2c 

VAD = • (HI) 

Thus, for typical conditions in a shielded core in which the ionization is dominated by cosmic rays, 
so that i/ff ~ 10, the ambipolar diffusion advances into the envelope at about 1/3 the speed of the 
expansion wave. The result is not changed significantly if we instead compare our result with that 
for the collapse of a thin disk with a frozen- in magnetic field: Li & Shu (1997) have shown that 
in this case the expansion wave advances at Cy/2, about the same velocity as in the unmagnetized 
case. 

The accretion rate is given by equation (101). Evaluating this rate at the time t' cv , we find 

M=£(^=5.4xW^r^r /2 m 9F -., ( 112 ) 



4 \t' AD J \5 M Q J V0.5 pc 

in good agreement with the estimate for the non-singular case in §5.7. If we express the accretion 
rate in terms of the sound speed required to support the sphere, as we did in equation ( p.ll| ), we 
find 

Recall that for the non-magnetized singular isothermal sphere, the accretion rate is M = 0.975c 3 /G 
(Shu 1977). Our results indicate that when ambipolar diffusion is included, the accretion rate 
is reduced only slightly: for z/g ~ 10, we obtain an accretion rate that is about 2/3 that for the 
equivalent non-magnetized case. 

This result is quite remarkable: it indicates that the presence of a strong magnetic field does 
not significantly reduce the accretion rate onto a protostar at late times in a cloud in which the 
ionization is due to cosmic rays (photoionization by FUV or X-ray photons increases and 
reduces the accretion rate proportionately). To see why this is so, let us first estimate the accretion 
rate for the singular isothermal sphere. Since about half the swept up mass has collapsed to the 
origin at any time, the accretion rate is Mgjg ~ M/2tg, where ts is the free fall time based on the 
mean density p = 3po. In terms of ts', the free fall time at po, we then have 



2 WJ 

This is somewhat higher than the exact answer, presumably because we neglected the retarding 
effects of the gas pressure. 

Next, consider the magnetized case. For a uniform cloud, the average accretion rate is 
M ~ M/2tAD based on the results in §4. Since the mass of a singular sphere is 3 times that of a 
uniform sphere with the same density at the surface, the ambipolar diffusion velocity, and hence 
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the collapse rate, will be three times greater: Msis.ad ~ \{M/txT>'), which is within a factor 1.5 
of the exact answer in equation (|112|). With these estimates, we find 



M SIS , AD _ 3M/2t AD > _ ^ / X = y/3 _ 



Msis V3M/2%' v V*ad7 / 8 \ 1/2 



With allowance for the fact that our estimate for this ratio is off by a factor 2 (the estimate for 
the magnetized case was somewhat too low, and that for the unmagnetized case was too high), 
this agrees with the result in equation ( |113| ). Several factors have raised the magnetized accretion 
rate close to the unmagnetized rate: in the magnetized case, the entire cloud is contracting at 
t = 2 £aDj whereas in the unmagnetized case the collapse begins only after the expansion wave 
reaches the point in question; thermal pressure retards the unmagnetized case somewhat (we 
have assumed that magnetic pressure dominates thermal pressure in the magnetized case); and 
finally, in going from an estimate based on a uniform sphere to a singular sphere, the rate for the 
magnetized case increases by p/po = 3, whereas that for the unmagnetized case increases only by 
the square root of this ratio. 

We conclude that, contrary to one's expectation (Mouschovias 1991), the accretion rate onto 
a protostar at late times after its formation is not significantly inhibited by the presence of a 
strong magnetic field, so long as the ionization is due to Galactic cosmic rays. Note that this 
result does not validate the singular isothermal model in toto. It only means that at late times, 
after the formation of the protostar, the mass accretion rate onto the central object approaches 
the limiting case of a singular isothermal cloud. We must stress, again, that our results for core 
collapse and immediately thereafter are at strong variance with the predictions of the singular 
isothermal sphere model. 

In our idealized model, the accretion will continue until the mass of the envelope is exhausted. 
In reality, several other factors could intervene: the outer envelope is likely to be photoionized 
(McKee 1989), with a much longer ambipolar diffusion time, and the outflow from the protostar 
could disperse the envelope (Shu et al 1987). 



6. Conclusions 

In this work we have investigated the evolution of a cold, magnetized molecular cloud by 
means of an idealized, semi-analytic model. We neglect thermal pressure, assume quasi-static 
evolution, and impose spherical symmetry. Our principal conclusions are the following: 

1. In the absence of thermal stresses, an initially homogeneous collapses homologously, without 
developing any structure. The collapse reaches a singularity at a time 2iAD> where £ad is 
the initial ambipolar diffusion time. 



-32- 



2. The central core of a centrally condensed cloud also collapses after a time 2t\D (where £ad 
is measured at the center of the cloud), leaving behind a quasi-static envelope. 

3. In a centrally condensed cloud the collapse of the core proceeds from outside-in, while the 
quasi-static envelope collapses from inside-out after the collapse of the central region. 

4. The pre-collapse evolution of the core is similar to the evolution of a homogeneous cloud. In 
particular, the quasi-static evolution of the central density is a universal function of time, 
p c oc (2 — t/tAD)~ 2 , independent of any parameters that enter the calculation. 

5. The core evolves quasi-statically up to the point that it becomes magnetically supercritical. 
Thereafter, the ions are no longer approximately static, and our approximations break down. 
As a result, we cannot treat the subsequent dynamical collapse. Nonetheless, even after the 
formation of the central protostar, we are able to follow the quasi-static evolution of the 
envelope so long as it remains magnetically subcritical. 

6. Most of the initial magnetic flux in the core is lost via ambipolar diffusion; after collapse, we 
estimate that the remaining trapped flux is about half the critical value. 

7. The mass accretion rate in the quasi-static envelope before and after core collapse is 
remarkably constant and of order 

*—» - 7m --' m 3 ' 2 m 3 ' 2 (mT' 2 ■ 

where Mo and Ro are, respectively, the cloud's mass and inital radius and e is the inital ratio 
of mean to central density, with 0.01^e^0.2. 

8. Much of the mass of the star is built up during the dynamical phase of accretion, which 
lasts about a free-fall time. Although we cannot treat this dynamical phase, our model does 
enable us to determine the total mass accreted and to follow the subsequent quasi-static 
accretion. Just as in the non- magnetic case (Foster & Chevalier 1993), the accretion rate 
during this stage is much larger than that of a singular isothermal sphere. 

9. The late-time mass accretion rate onto the protostar is quasi-static. Notwithstanding the 
inhibiting effects of the magnetic field, the accretion rate becomes comparable to that of a 
singular isothermal sphere. 

10. The excellent agreement between our results for the pre-collapse evolution and the detailed 
numerical work of Fiedler & Mouschovias (1992) indicates that our simple approach captures 
the basic physics relevant to the evolution of a cloud undergoing ambipolar diffusion. Our 
results are especially relevant to the late stages of cloud evolution, when thermal stresses 
are relatively small. Since thermal pressure does play an important role at earlier times, 
particularly in creating the initial density distribution, our intention is to extend the present 
calculation to one that incorporates thermal pressure in a subsequent paper. 
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Fig. 1. — Evolution of the Lagrangian radius r (in units of the initial cloud radius i?o) as a function 
of time for a homogeneous cloud with 5 = 4 x 10" 3 . The five curves show the evolution of r for 
five different values of the dimensionless Lagrangian coordinate m: m = 10 -2 , 3.16 x 10~ 2 , 10 _1 , 
3.16 x 1CP 1 , and 1. All the masses are in units of the cloud's mass Mq, and time is in units of the 
characteristic ambipolar diffusion time, £ad- 

Fig. 2. — Evoution of the magnetic critical ratio M$/M for a homogeneous cloud with 5 = 4 x 10 -3 . 
(a): M$/M as a function of time (in units of the ambipolar diffusion time £ad) for three different 
values of the dimensionless Lagrangian coordinate m: m = 10~ 2 , m = 10 -1 , and m = 1.0. (b): 
M<$>/M as a function of the dimensionless coordinate M/Mq at three different times (in units of 
the characteristic ambipolar diffusion time): r = 0.0, 1.62, and 1.86. 

Fig. 3. — Ratio of the ion's velocity Vi to the neutral's velocity v n for a homogeneous cloud with 
5 = 4 x 10 -3 . Shown is Vi/v n as a function of the dimensionless Lagrangian coordinate M/Mq at 
four different times (in units of the characteristic ambipolar diffusion time): r = 0.0, 1.00, 1.62, 
and 1.86. 

Fig. 4. — The mass of the initial central core, M c (in units of Mo), as a function of the initial degree 
of concentration, e = (p) / p c $ for a stratified cloud with p = 2. Only values m c = M c /Mq < 1 are 
physically relevant. Clouds with e > 3/7 have m c > 1, and therefore collapse as a homogeneous 
cloud. 

Fig. 5. — Lagrangian radius r, neutral velocity v n , density p, and magnetic field B for our fiducial 
stratified cloud with an initial degree of concentration e = 0.1 and magnetic parameter 5 = 4 x 10~ 3 
at five different times (in units of the characteristic ambipolar diffusion time £ad) : t = 0.0, 0.88, 
1.37, 1.64, and 1.80 (between consecutive values of r the central density increases by a factor of 
(10) 1 / 2 ). (a): Lagrangian radius r (in units of the cloud's initial radius Ro-top panel) and neutral 
velocity v n (in units of the fiducial velocity v n fl~bottom panel) as a function of the dimensionless 
Lagrangian coordinate M/Mq. Time increases from top to bottom for the r-curves and from bottom 
to top for v n . The dashed line indicates the location of the central core, (b): Density in units 
of the initial central density p Ct o (top panel) and magnetic field (plotted as B/Bq — 1, where Bq 
is the strength of the magnetic field at the edge of the cloud; bottom panel) as functions of the 
dimensionless Lagrangian coordinate M/Mq. Time increases from bottom to top. The dashed line 
indicates the location of the central core. 

Fig. 6. — The acceleration of the neutrals in units of the local gravitational acceleration. Shown 
is a, for our fiducial cloud, as a function of the dimensionless Lagrangian coordinate M/Mq and 
dimensionless reversed time (measured from the formation of a singularity at the origin) 2 — t/t ad ■ 

Fig. 7. — Time evolution of the magnetic flux $ in our fiducial cloud (in units of the fiducial 
magnetic flux ttBq Rq) as a function of the dimensionless Lagrangian coordinate M/Mq (left panel) 
and the Eulerian radial coordinate R (in units of the cloud's initial radius Rq] right panel). Each 
curve corresponds to a different dimensionless time r ( r = 0.0, 0.88, 1.37, 1.64, and 1.80), which 
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increases from top to bottom in the left panel and from bottom to top in the right panel. The dashed 
line and the arrows indicate the location of the central core; the arrows are labeled by the value of 
t/t ad- 

Fig. 8. — Ratio of the ion to neutral velocity, Vi/v n for our fiducial cloud (e = 0.1 and S = 4 x 10~ 3 ), 
as a function of the dimensionless Lagrangian coordinate M/Mq for five different times (in units of 
the characteristic ambipolar diffusion time £ad; from bottom to top): r = 0.0, 0.88, 1.37, 1.64, and 
1.80. The dashed line indicates the location of the central core 

Fig. 9. — Evolution of the critical ratio M<j>/M for our fiducial cloud (e = 0.1 and 5 = 4 x 10~ 3 ). 
M$/M is plotted as a function of the dimensionless Lagrangian coordinate M/Mq for five different 
times (in units of the characteristic ambipolar diffusion time £ad! from top to bottom): r = 0.0, 
0.88, 1.37, 1.64, and 1.80. 

Fig. 10. — Ratio of thermal to magnetic stresses, /?', as a function of the dimensionless Lagrangian 
coordinate M/Mq and dimensionless reversed time (measured from the formation of a singularity 
at the origin) 2 — t/t ad- Shown are the results for our fiducial cloud, e = 0.1 and 5 = 4 x 10~ 3 , 
with initial central number density n Cj o = 2200 cm~ 3 , magnetic field strength Bq = 30 /xG at the 
edge of the cloud, and temperature T = 10K. 

Fig. 11. — Comparison with the numerical results of Fiedler & Mouschovias (1992) for the time 
evolution of the central number density n c (top panel) and central magnetic field B c (bottom panel) 
in a cloud with mass Mq = 45.5 M Q , initial radius Rq = 0.75 pc, and external magnetic field 
Bq = 30/^G. Shown are the results of Fiedler &; Mouschovias for their Model 1 (filled circles and 
dashed line) after the initial relaxation and up to the onset of dynamical contraction and those 
from our semi-analytic model (solid line). 

Fig. 12. — Comparison with the numerical results of Fiedler & Mouschovias (1992) for the evolution 
of the density profile. Shown is the number density as a function of dimensionless radius v/Rq for 
three different times during the quasi-static phase. The filled circles are the numerical results for the 
midplane of Fiedler & Mouschovias's Model 1 (light grey: t = 10.2 Myr, dark grey: t = 15.1 Myr, 
and black: t = 16.0 Myr); the solid lines are the results from our semi-analytic model. For 
easier visualization, the curves have been shifted horizontally by a constant C: C = log 0.5 for 
t = 10.0 Myr; C = log 0.75 for t = 15.1 Myr; and C = 0.0 for t = 16.0 Myr. 

Fig. 13. — Core's magnetic flux $ C orc (in units of the cloud's initial magnetic flux $o) as a function 
of dimensionless reversed time (measured from the formation of a singularity at the origin) 2—t/t\j). 
Shown is $ core /$o for 5 = 10 -3 and four different values of e: e = 10~ 4 , 10~ 3 , 10~ 2 , and 10 _1 . 

Fig. 14. — Post-collapse evolution of the density profile (top panel) and neutral velocity profile 
(bottom panel) for our fiducial cloud (e = 0.1). Shown are the density in units of the initial central 
density and the neutral velocity in units of the fiducial velocity v Ht o as a function of radius in the 
cloud (in units of the initial cloud radius Rq) for six different times: r = 2.0, 2.1, 2.2, 2.3, 2.4, and 
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2.5. Time increases from the topmost curve down. The inner boundary of each curve corresponds 
to the point where the acceleration of the neutrals is 10% of the local gravitational acceleration. 

Fig. 15. — Mass accretion rates as a function of radius for our fiducial cloud (e = 0.1). Shown is 
M (in units of the fiducial mass accretion rate Mad) as a function of radius in the cloud (in units 
of i?o) before the collapse of the inner core (dotted lines) and afterwards (solid line and shaded 
region). For the pre-collapse accretion rates each curve corresponds to a different time (in units 
of the characteristic ambipolar diffusion time £ad; from bottom to top): r = 0.0, 0.88, 1.37, 1.64, 
and 1.80. The solid line is the accretion rate at the time when the core collapses (r = 2.0) and 
the shaded region corresponds to the post-collapse mass accretion rates for r < 2.5; the individual 
curves in this regime are too close to each other to show clearly on the scale of this plot. The 
post-collapse curves (r > 2.0) are truncated to the left at the point where the acceleration of the 
neutrals is 10% of the local gravitational acceleration. 

Fig. 16. — Protostellar mass M* and mass accretion rate onto the protostar, M*, as functions of 
dimensionless time (measured from the formation of a singularity at the origin) t/t^D — 2. Shown 
are M* (in units of the cloud's mass Mo; top panel) and M* (in units of the fiducial mass accretion 
rate Mad; bottom panel) for three different values of e: e = 0.05, 0.1, and 0.2. 
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